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A striking feature of the marine ecosystem is the regularity in its size spectrum: the abundance 
of organisms as a function of their weight approximately follows a power law over almost ten orders 
of magnitude. We interpret this as evidence that the population dynamics in the ocean is approxi- 
mately scale-invariant. We use this invariance in the construction and solution of a size-structured 
dynamical population model. 

Starting from a Markov model encoding the basic processes of predation, reproduction, mainte- 
nance respiration and intrinsic mortality, we derive a partial integro-differential equation describing 
the dependence of abundance on weight and time. Our model represents an extension of the jump- 
growth model and hence also of earlier models based on the McKendrick-von Foerster equation. 
The model is scale-invariant provided the rate functions of the stochastic processes have certain 
scaling properties. 

We determine the steady-state power law solution, whose exponent is determined by the relative 
scaling between the rates of the density-dependent processes (predation) and the rates of the density- 
independent processes (reproduction, maintenance, mortality). We study the stability of the steady- 
state against small perturbations and find that inclusion of maintenance respiration and reproduction 
in the model has a strong stabilising effect. Furthermore, the steady state is unstable against a 
change in the overall population density unless the reproduction rate exceeds a certain threshold. 



I. INTRODUCTION 

The population dynamics in the pelagic marine ecosys- 
tem (the open sea) are particularly amenable to mathe- 
matical modelling and analytic understanding. Because 
size is the most important factor in determining who eats 
who, rather than species pQ, it is possible to work with 
a model in which a single function </>(«;, t) describes the 
total population density of organisms of weight w at time 
t, aggregated over all species. We can build on many pre- 
vious works using partial integro-differential equations to 
describe the time evolution of this total population den- 
sity, in particular [2HS]. 

In this paper we exploit an additional special property 
of the pelagic zone: its approximate physical scale in- 
variance. The open sea looks similar at a range of scales. 
That is due to the fact that over many orders of magni- 
tude there are no physical features and no strong phys- 
ical principles that would single out a particular inter- 
mediate scale. We expect that this approximate scale 
invariance of the environment breaks down only at small 
scales at which molecular diffusion processes dominate 
and at large scale at which geography affects ocean cur- 
rents. This kind of scale-invariance is not present in ter- 
restrial environments, where not only the physical struc- 
ture sets a scale, but where in addition the effects of 
gravity quickly become important for larger organisms. 
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It is not a priori guaranteed that the scale invariance 
of the physical environment will also lead to scale invari- 
ance of the ecosystem. The organisms populating the 
environment and their interactions could break the scale 
invariance. This paper is however based on the assump- 
tion that evolution, in its drive to make use of all available 
ecological niches, has made optimal use of the scale in- 
variant environment by filling it with an ecosystem that 
roughly preserves scale invariance. It is beyond the scope 
of the current paper to investigate the evolutionary mech- 
anisms that would lead to the self-organisation of such a 
scale invariant ecosystem. We will be content with deriv- 
ing the consequences of scale invariance for the dynamics 
of the population density. 

The main observational evidence for approximate scale 
invariance is that the equilibrium size distribution of or- 
ganisms in the open ocean is approximately given by a 
power law 4>{w) oc w -7 , valid over almost ten orders of 
magnitude [10Hl5] . Several theoretical models have been 
proposed in the past to derive this steady-state size spec- 
trum [21 131 [MSI [23 HZ] but, as far as we know, in this 
paper we are the first to point out that this can be un- 
derstood as a consequence of the scale invariance of the 
underlying dynamical model. 

The main processes that affect the abundance of organ- 
isms as a function of weight are predation, reproduction, 
maintenance respiration and intrinsic death. We will in 
the next section write down a dynamical model that in- 
corporates all these processes. In general such a model 
will not reproduce the observed power-law size-spectrum 
in the steady state. This is the reason why existing work 
either only models predation [21 [31 [SHE] or makes the sim- 
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plifying assumption that reproduction, maintenance res- 
piration and death are exactly proportional to predation 
[5]. We are able to avoid those difficulties by exploiting 
scale invariance. 

While scale invariance leads to a power-law steady- 
state solution, it docs not guarantee that this steady 
state is a stable equilibrium solution. Indeed, a linear 
stability analysis of the jump-growth model [18 , which 
is also scale invariant, showed that, for realistic choices 
of the parameters, the steady state in that model is un- 
stable against small perturbations. The steady state ob- 
served in nature is stable, so the pure jump-growth model 
is missing some important stabilising effect. This was 
the motivation for the investigations of the more general 
model in this paper. We will see that the inclusion of 
maintenance respiration and reproduction has a strong 
stabilising effect. 

We do not model spatial variation in the population 
density. Nevertheless our model formally resembles a 
model in one space and one time dimension, where the 
space dimension represents the weight of the organisms. 
This formal analogy between space and weight will be- 
come particularly clear when we transform from the 
weight variable w to the logarithmic weight variable x 
in Section HVl A scale transformation in w then corre- 
sponds to a translation in i. In the special case where 
the scale transformations do not affect the time variable 
we end up with a model that is translationally invariant 
in both space and time, and many of the usual techniques 
can be applied, like for example the use of the standard 
Fourier transform in the solution of the linearised model 
in Section |VH 

Another difference between our model and those usu- 
ally studied by physicists derives from the fact that feed- 
ing is a non-local interaction in weight space. Most fish 
and plankton do not feed on other individuals that are 
close to their own weight. Instead they prefer prey that 
are substantially smaller. Similarly they produce off- 
spring at a weight far below their own. This is the 
reason why the population density <p(w,t) is modelled 
not by a partial differential equation but by the partial 
integro-differential equation ([8]) involving integral terms 
that encode the feeding behaviour (|9| and reproductive 
behaviour ( |10[ ). In the linear stability analysis this leads 
to the non-local dispersion relation ( |49| . 

The organisation of this paper is as follows. In Sec- 
tion [Tl] we formulate a stochastic model encoding feeding, 
reproduction, maintenance, and death. We do not make 
any assumptions about the rates for these processes but 
keep them as general parameters of the model. By taking 
the macroscopic limit of the stochastic model we derive 
the deterministic evolution equation ([8]). The derivation 
of a Fokker-Planck equation for the stochastic fluctua- 
tions away from the macroscopic model is left to Ap- 



use scale invariance to find a power-law steady state so- 
lution of the model in Section [Vj where we discuss as well 
its main properties and restrictions. The overall popula- 
tion level is predicted by the model and we show in Ap- 
pendix [B] that is positive under biologically relevant as- 
sumptions. A time-damped power-law solution can also 
be found, as we show in Appendix [C] In order to in- 
vestigate the stability of the steady state we determine 
the spectrum of small perturbations around the steady 
state in Section IVll We summarise our results in the Sec- 
tion (VTTJ In Appendix [D] we describe how variability in 
the absorption efficiency can be approximated by the ad- 
dition of a diffusion term to the model. 



II. SIZE-STRUCTURED POPULATION MODEL 

We start by constructing a stochastic model for the 
dependence of abundance of individuals on weight and 
time, taking into account the basic processes of preda- 
tion, reproduction, maintenance and intrinsic mortality. 

As in previous work [5] , instead of keeping track of the 
weight of each individual, we aggregate individuals of 
similar weight into discretised weight brackets [wi,Wi+i) 
for i e Z. Weight is the only attribute of individuals 
that is used in this model. Species identity and life 
stage are ignored. The weight distribution of individ- 
uals in a large fixed volume Q is described by the vector 
n = (. . . , no, 71%, . . . ) whose entries give the number 
of organisms in each weight bracket. We will later let the 
size of these brackets go to zero to obtain the continuum 
model. Working with discrete weight brackets allows us 
to avoid the mathematical complexities involved when 
trying to describe a stochastic model directly in contin- 
uous space, see for example |19J . 

The primary stochastic processes involved in the model 
are illustrated in Figure [T] which shows the various ways 
in which an event can affect the number of organisms 
in bracket i. Below we will show how the deterministic 
equation can be read off directly from this figure. In Ap- 
pendix [S] we start from a master equation for the prob- 
ability distribution P(n,t) and carry out the systematic 
expansion of van Kampen |20j in powers of the inverse 
system volume f2 _1 . The same method was used in [5] to 
derive a jump-growth equation. This method is based on 
splitting each variable rii(t) into a deterministic, macro- 
scopic component 4>i (t) describing the density of individ- 
uals in weight bracket i, and a fluctuation component 

m(t) as 



(1) 



pendixjAj In Section III we impose scale-invariance and 
derive sufficient scaling conditions (15 1 for the rate func- 



tions. Section[TV]is devoted to transforming the evolution 



equation into the more convenient form (27). We then 



The powers of volume are chosen so that the new vari- 
ables <f>i and rji no longer scale with the system volume. 
This method not only gives the macroscopic behaviour 
for the densities <f>i(t) at leading order in the expansion, 
at higher orders it describes the stochastic fluctuations 
around the macroscopic solution as well, giving at next- 
to-leading order a linear Fokker-Planck equation. How- 
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FIG. 1. Individual stochastic processes involved in the model. 
The figure shows all the events that affect the number of or- 
ganisms in weight bracket i. The arrows indicate the move- 
ment of individuals between weight brackets, (a) Predation: 
There are three predation processes affecting bracket i. The 
first represents a predator in i eating a prey in j and absorb- 
ing enough prey weight to end up in k. The second represents 
a prey in i being eaten and the third represents a predator 
entering i after feeding, (b) Reproduction: in the first pro- 
cess an individual in i produces rriijk offspring in bracket k 
and, as a consequence, decreases its weight to j. In addition, 
there are two processes that increase the number of individ- 
uals in i. (c) Maintenance respiration: these processes move 
individuals to the next-lower bracket. There are two ways in 
which this affects to bracket i. (d) Intrinsic mortality: with a 
certain rate, a single individual in bracket i is removed from 
the system. 



ever, because the system volume is so large, in this pa- 
per we will concentrate on the macroscopic, deterministic 
equation. 

From Figure [l] we can obtain the contributions to the 
time evolution of fa from each of the processes we con- 
sider: predation (P), reproduction (B for Birth), mainte- 
nance respiration (R for Respiration) and intrinsic mor- 
tality (D for Death), 



dt 



dt 



dt 



dt 



dfa 
"dt 



(2) 



P V / B V / R 

We will describe each of them in turn below. 

A. Predation 

A predation event moves a predator from a weight 
bracket i before feeding to a higher weight bracket k af- 
ter feeding and removes a prey organism from a weight 
bracket j. Let Pijk be the rate constant for such a pre- 
dation event. As illustrated in Figure [Ib, there are three 
ways in which a predation event can affect the number 
7i{ of individuals in a weight bracket i: 1) an organism 



in bracket i can eat another organism and grow; 2) an 
organism in bracket i can be eaten by another organism; 
and 3) an individual in a lower bracket can absorb enough 
prey weight and grow into bracket i. Because the rate at 
which a particular individual will encounter prey will be 
proportional to the density of prey, the probability that 
one of the Uj individuals in bracket i eats any of the rij 
individuals of bracket j and increases its size to bracket 
k in the time interval from t to t + dt is Piji~ninjVL~ 1 dt. 
Hence the contribution of the predation events to the de- 
terministic time evolution of the density fa of organisms 
in bracket i is 



dt 



Pijkfafa - Pjikfafa + Pjkifafa)- (3) 



Our model for predation can be viewed as a generali- 
sation of the Smoluchowski coagulation equation [3TJ [55] 
which is obtained in the special case when the resulting 
weight Wk is equal to the sum of the weights Wi and Wj 
of predator and prey. 



B. Reproduction 

Most fish reproduce by laying a large number of eggs 
that are subject to heavy predation. Only a small frac- 
tion of the eggs survives to hatch and join the consumer 
spectrum. In a size-spectrum model however, in which 
size is the only attribute of an organism and its life stage 
is ignored, all organisms are assumed to be prey and 
predator simultaneously from the moment they are born. 
We can therefore not model the egg life stage and there- 
fore can not provide an entirely realistic model of the 
reproductive processes. Our model is thus designed to 
only capture two features of reproduction that we deem 
essential to the size-spectrum dynamics, namely that it 
moves biomass up the size spectrum from large weight 
to small weight and that it replenishes the population 
numbers at smaller weight. 

We assume that a reproduction event moves a par- 
ent organism from a weight bracket i to a lower weight 
bracket j and produces a large number rriij}. of smaller 
offspring in weight bracket k. We set the number of off- 
spring to niijk — [(wi — Wj)/wk\, where [^J denotes 
the largest integer smaller than x, so that their com- 
bined weight is approximately equal to the weight lost 
by the parent. Let Bijk be the rate constant for such 
an event. The probability that one of the rii individuals 
in bracket i reproduces in such an event in dt is then 
Bij^riidt (note that reproduction, unlike predation, is a 
density-independent process) . 

As depicted in Figure[T]D, there are three ways in which 
reproduction changes the number of the organisms in 
bracket i: 1) the parent belongs to the i-th bracket, and 
after spawning moves to a lower bracket; 2) a parent 
looses weight during spawning and moves to bracket i; 
and 3) the bracket i receives the offspring from a repro- 
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duction event. The contribution of these three possibili- 
ties to the deterministic equation is therefore given by 



dt 



^ {-Bijk<pi + Bj ik (f>j + rrijkiBjki^j) . (4) 



C. Maintenance and mortality 

In between feeding events, organisms continuously 
draw upon their reserves to maintain themselves. The 
weight loss due to maintenance respiration is modelled 
by events that move individuals to the next-lower weight 
bracket, assuming that the width of each interval is small 
enough (this is not a restriction, because at the end we 
will take the continuum limit where all these widths tend 
to zero). The probability that any of the n, individ- 
uals in the bracket i undergoes such an event in dt is 
RiUidt, where Ri is the maintenance respiration rate for 
the bracket i. Figure [TJ; shows the two ways that these 
primary processes change the number of individuals in 
bracket i, giving the contribution 



dt 



= R 



i+l<Pi+ 



i — Ri 



(5) 



Organisms can also die for reasons other than being 
eaten. This introduces a fourth process in the model 
that accounts for intrinsic mortality. With a probability 
DiUidt, a single individual in bracket i is removed from 



the system in dt (see Figure [TJi) . This gives the contri- 
bution 



dt 



-D, 



(G) 



Note that, like reproduction, both maintenance and in- 
trinsic mortality are modelled as density-independent 
processes. 



D. Continuum limit 

We now take the continuum limit of each contribution 
to the macroscopic equation by writing A, = Wi + \ — Wi 
and then taking the limit Aj — > uniformly for all i. 
The variables </>i(i) are combined into a density 4>{w, t) of 
individuals per unit weight per unit volume as a function 
of weight and time so that 4>(wi, t) = 4>i{t)/ A^. The sum 
over weight brackets is replaced by an integral, J^. Aj — > 
J dw. Continuum rate functions are introduced as 



P(w t ,W jl W k ) = Py fe /A fe , 

B(wi,Wj,w k ) = B ijk /(AjA k ), 
R{wi) = AA, 
D{w l )=D l . 



(7) 



After taking the continuum limit, the ordinary differen- 
tial equations ^ for the & assemble into a partial dif- 
ferential equation for (j)(w,t), 



d<Kw,t) = ( d<Kw,t) \ ( d</>(w,t) 
at \ dt J p \ dt 



d<f)(w,t)^\ ( f d(j)(w,t) 

R 



dt 



V dt 



(8) 



which will be the object of study in the remainder of the paper. The contribution from predation ^ in this limit 
becomes 

Integrals run over all positive weights. Similarly, the continuum limit of Q is 

' d<j)(w, t) 



dt 



= I dw' I dw" [-B(w,w',w")(t){w,t)+B(w',w,w")(j){w',t) + - —B{w',w",w)(j)(w',t)). (10) 



Maintenance and intrinsic mortality are described by 
/ d<j>(w,t)\ d 



V dt 



dw 



[R(w)<t>(w,t)] 



and 



d(f)(w,t) 
dt 



= -D(w)<t>(w,t), 



respectively. 



(11) 



(12) 



We note that the model with only the predation 
term reduces to the jump-growth equation derived in 
[5] when the feeding rate is chosen as P(w,w' ,w") — 
k(w,w')S(w + Kw 1 — w"), for some feeding preference 
function k(w,w') and a fraction K accounting for the 
feeding efficiency. As pointed out in [5], when the typi- 
cal predatonprey mass ratio is sufficiently large and the 
population density is sufficiently smooth that model in 
turn can be approximated by the McKendrick-von Foer- 
ster equation [23 , 24J . It is this simplified model in terms 
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of the McKendrick-von Foerster equation that forms the 
basis of most previous analytic studies of the size distri- 
bution [2H7]. We will see that, thanks to scale invariance, 
the much more general evolution equation ^ is amenable 
to similar analytical investigations. 



III. SCALE INVARIANCE 

We now derive restrictions on the feeding rate P, the 
reproductive rate B, the maintenance rate R and the 
mortality rate D that guarantee the scale invariance of 
our model. 

As we have a model describing the dependence of the 
abundance on two variables, namely weight and time, 
we can consider separate scaling in weight and in time. 
However, we expect our model to be invariant only under 
simultaneous scaling of both weight and time, because we 
expect life processes to run faster for smaller organisms. 
Hence, when we change the weight scale by a factor c we 
should also change the time scale by a factor (fi, where the 
constant £ expresses how the speed of the dynamics scales 
with weight. So we will consider the transformations 



(w, t) h-> (cw, c^t) with c > 0. 



(13) 



Under such scale transformations the density <p(w, t) 
transforms as 



(w, t) h-> c 7 0(ciu, c£t) 



(14) 



with some, so far undetermined, exponent 7, also called 
the scaling dimension of <j> [25] . 

Requiring the evolution equation Q to remain un- 
changed under this transformation imposes the condi- 
tions 



c~^~ 2 B(w,w' ,w") = B(cw,cw' ,cw"), 
c^Riw) = R[cw), 
c~^D(w) = D(cw). 

We choose to factorise the feeding rate as 

P(w, w', w") = S(w, w')A(w, it/, to"), 



(15) 



(16) 



where S(w, w') determines the rate at which a predator 
of weigh w eats a prey of weight w' and A(w, to', to") gives 
the probability density that such a feeding event makes 
the predator grow to weight to". On average only a cer- 
tain proportion Q of the prey biomass will be absorbed by 
the predator, so that on average w" = w + Qw', but there 
will be some variability, due to differences in predator di- 
gestion and prey composition. This will be modelled by 
the probability density A(w, to', to"), so that 



dw" w"A(w, u/, to") = w + Qw'. 



(17) 



We should make it clear that the model, and therefore 
the results, depend only on the combined rate P. We 
factorise it into S and A only to make the ecological 
origin of the rate clearer [5J [8] , but how we choose this 
factorisation has no influence on the results. 

Because the model only depends on the product of S 
and A, we are free to choose the relative scaling between 
these factors, as long as the product scales as in ( 15 ). We 



want A to scale as a probability density in to" and hence 
we choose 



A(w,w',w") = w"- 1 A (w/w',w"/w), 



S(w, to') = (w/w ) 7 1 S (w/w'), 



(18) 



Here wq is an arbitrarily chosen reference weight. We 
have also introduced the scaling functions Aq and So that 
are invariant under scale transformations. 

Similarly, we require that the reproduction rate B scale 
as a density in to' and to", so according to the scaling rules 
given in (151 we can write 



B(w, w', w") = w'- 1 ™"- 1 (w/w Q ) e B (w/w', w"/w) 

(19) 

for some scaling function B . This behaviour for the 
reproduction rate implies that the average weight of an 
offspring is proportional to to 1 ^, where w is the weight 
of its parent. This should however not be taken as a 
prediction of the scaling of egg sizes because, as stressed 
previously, our model of reproduction is not formulated 
at a sufficient level of detail for that purpose. 

The scale transformations for the maintenance and 



death rates given in (151 allow us to express them as 



R(w) = (w/wq) 1 i R(w ), 
D(w) = {w/wo)^D(w ), 



(20) 



The functions So, Aq and Bq, the constants R(wq) and 
D(wq) and the exponents 7 and £ are not fixed by the re- 
quirement of scale invariance and need to be determined 
from observations or separate theoretical arguments. 

The restrictions that we impose to achieve scale in- 
variance predict the relative scaling of these rates. In 
particular, the scaling of all the rates contains the ex- 
ponent £. It is widely believed that the maintenance 
metabolic rate scales as w 3 ^ 4 [2"6T - [2"8] . Comparing that to 



the scaling in (20 1 gives £ = 1/4. This then predicts that 
the mortality rate scales as w^ 1 / 4 , which is in line with 
observations [28l[29]. Moreover, the exponent 7 of the 
scaling transformation ( 14 ) of the density is determined 



once the scaling of the feeding rate and £ are known. 
If we assume, following Benoit and Rochet [BJ, that the 
feeding rate is proportional to the volume searched per 
unit time, we can identify the scaling exponent of S as 
the exponent of the search volume, which is around 0.8 
[30] . Assuming £ = 1/4 it follows that 7 « 2, which 
is in line with observations as well [9]. We will see in 
Section [V] that the exponent 7 coincides with the expo- 
nent of the power-law steady-state solution. Therefore, 
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the relative scaling between density-dependent processes 
(feeding) and density-independent processes (reproduc- 
tion, maintenance and death) determines the steady state 
exponent. 

Besides scale invariance, our model also has time trans- 
lation invariance. Time translations and scale trans- 
formations together generate the non-abelian symmetry 
group Aff(l,R), the group of orientation-preserving lin- 
ear transformations of the real line, also known as the 
ax + b group. In the special case £ = the symme- 
try group becomes abelian, a point which we will ex- 
ploit in Section \VJ\ because in that case the equation is 
translationally invariant and hence the standard Fourier 
transform can be applied to solve the equation for small 
perturbations of the steady-state. 

The symmetry is enhanced when the model only con- 
tains the predation terms, which are quadratic in (f>. In 
that case there is also an invariance under scaling in time 
alone, 



(w,t) H> \4>(w,\t) 



(21) 



IV. CHANGE TO LOGARITHMIC WEIGHT 

For the upcoming analysis it is convenient to make a 
change of variable to x = log(w/wo), where wq is the 
arbitrarily chosen reference weight. We refer to x as the 
logweight. A scale transformation w i— ► cw corresponds 
to a translation x <— > x + log c. 

We introduce the density u{x) so that Qu(x)dx is the 
number of individuals in volume ft with a logweight be- 
tween x and x + dx. Thus u(x) — w<p(w) and we can 
easily translate our results for u(x) back into results for 
<j)(w), if desired. Under a scale transformation the func- 
tion u(x,t) transforms to c 7 ~ 1 u(a; + log c,(fit). We will 
apply this change to the various terms of the evolution 
equation (si). 

We first introduce the predation rate p such that 
p(x,x' ,x") — w" P(w,w' ,w"). The factorisation of P 
into S and A in equation ( 16 1 together with equation 
(181) leads to 



for any A £ 



p(x, x' , x") = e^ p ^' x s(x — x')a(x — x' , x" — x), (22) 

where we have defined the functions s{y) = So(e v ) and 
a(y, z) — Ao(e v , e z ) and we have introduced the exponent 
p = 7 — 1 for latter convenience. We now transform to 
logweights, substitute this form for p into (9| and perform 
a change of variables in each of the feeding terms so that 
the integration variable coincides with the argument of 
s. The result is 



' du(x) 
dt 



(23) 

where we have taken into account the fact that a(y, z) is a probability density and hence it is normalised. The integrals 
all run over the whole real line. In what follows, we will often not indicate the time-dependence of u(x,t) explicitly 
but write just u(x) instead, as in the above equation. 



We can perform the same changes for the reproduction term. Using the scaling form (19) for B in ([10]), then 
transforming to logweights and defining b(y, z) — Bo(e y ,e z ), we obtain 



^^M^j =e~ ix j dy J dz b{y, z)( - u{x) + e~ iy u(x + y) + e ( «~ 1)z (l - e - y )u(x - z)). 



(24) 



Finally, the contribution of maintenance and intrinsic 
mortality can be expressed in logarithmic weights as 

(^) < 25 » 



and 



du(x) 
dt 



= -de ix u(x), 



(26) 



respectively. Here we have introduced the constants r 
R(wq)/wo and d = D(wq). 



The full dynamical equation in terms of logweights is 



du 
dt 



du 



du 
~dt 



du 

~dt 



du 



(27) 



We have chosen not to fully non-dimensionalise this evo- 
lution equation: t still has dimension of time, u(x) has 
the dimension of inverse volume, r, d and b have dimen- 
sion of inverse time, s has dimension of volume over time 
and a is dimensionless. 
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V. POWER-LAW STEADY STATE SOLUTION 



Conservation of number of individuals 



Solving the integro-diffcrcntial evolution equation (27) 



is difficult in general. However we can simplify the task 
by looking for solutions that are invariant under symme- 
try transformations. 

In this section we will study the solution that is in- 
variant under both time-translations and scale trans- 
formations. We leave the study of the general scale- 
invariant solutions to Appendix [C] Invariance under 
time-translations means that we arc looking for a steady- 
state solution <p{x) that has no dependence on time. In- 
variance under scale transformations ( 13 ) and ( 14 1 then 
implies 



(w) = 4>(w ) 



(28) 



where wq is the arbitrarily chosen reference weight. After 
transforming to logweights as in Section |TV| this takes the 
form 



u(x) = u^e" 



px 



(29) 



where uq — wq4>(wo) and p = 7 — 1. Substituting this 
form for the solution into the evolution equation (27) 
gives an equation for the overall population level Uq, 

c p uq = -r(p + £) - d + c B , (30) 

where the constants cp and cb are given by 

* = /*.M(e- + ^-e»/«,e<^ to , 2 )) 



CB 



dy / dz b(y,z)(-l + 



(1 



When cp ^ this uniquely determines uq and hence 
the steady state solution. We will show in Appendix [B] 
that uq is positive under some ecologically reasonable 
assumptions about the parameter functions. 

Note that scale invariance fixes the power-law form of 
the steady-state size spectrum and the steady-state ex- 
ponent p is determined entirely by the scaling behaviour 
of the parameter functions, see (15). It is not dependent 



on any other details of the interactions in the model. 

A special situation arises in the case where mainte- 
nance, reproduction and intrinsic mortality are absent 
from the model. In this case only the first scaling rela- 



tion in ( 15 ) remains and it is not enough to determine 
both £ and p. However an equation for p is obtained 



by noticing that in this case the right-hand side of (30) 
is zero and for uq ^ this implies that cp = 0. This 
constraint should then be used to determine the scaling 
exponent p given a particular choice for s{y) and a(y, z). 
The overall population level uq is not determined by the 
model in this case. This special situation was considered 
in most previous work 2 . 3 , 6 , 8 . 



There is a continual flux of individuals from lower 
weight to larger weight to make up for the losses due to 
predation and intrinsic death. In previous models that 
considered only the predation process [51 [5] there was no 
source for this influx of small individuals. Instead they 
appeared from x = —00. Now that we are modelling the 
reproduction process, we do have a source of individu- 
als and can impose that in the steady state this source 
should exactly balance the losses. 

The easiest way to impose this balance is to impose for 
each weight bracket i that the number of individuals en- 
tering the bracket from the left due to predation exactly 
equals the number of individuals leaving that bracket to 
the left, either as offspring or through weight-loss. This 
gives 

^ Pjki<t>j<t>k = /Xmjjk + ~t)B ijk <j)i + Ri4>i. (32) 

j,k j,k 

Thanks to scale invariance, all these conditions for differ- 
ent i are equivalent. In the continuum, after substituting 
the steady state solution and changing to logweight no- 
tation, this condition reads 



fpuo = /b + r, 
where we have defined the constants 

f P = [dy [dzs(y)a(y,z)e^ z+ Py, 



(33) 



(34) 



We can use this constraint, together with the steady-state 
condition ( 30 ) to fix the maintenance rate in terms of the 



/ B = dy dz b(y, z) (e- z (l - e"») + l) 



(31) other parameters of the model, 



(cb -d)f P - c P / E 
cp + (p + Ofp 



(35) 



If we use this to eliminate r from the equation (30) for 
Mo we get 



u 



c B -d+{p + Q/b 
cp + (p + 0/p 



(36) 



Obviously, for the model to make sense we need r to 
be positive. This is not possible for all choices of the 
other parameters. In particular, this requirement defines 
a range of allowed exponents p. To investigate this fur- 
ther, we will now make specific choices for the functions 
that appear in the feeding and reproduction rates. 



B. Choice of parameter functions 

For the prey selection function s(y) we choose a Gaus- 
sian that expresses that there is a preferred value j3 for 



the log of the predatonprey mass ratio and a certain vari- 
ance a a around this mean 1311. So we set 



s(y) = s ga Ay- P) 



with 



-x 2 /2<7 2 



2ircr 



(37) 



(38) 



The parameter sq has dimension of volume over time and 
sets the overall feeding rate. 

For the absorption probability density a the simplest 
assumption would be that a fixed proportion Q of prey 
mass is absorbed in all feeding events, i.e., that in terms 
of the predator mass w and the prey mass w' the mass 
after feeding is always w" = w + Qw' . This corresponds 
to a choice a{y, z) = 5(z — tp(y)) where 



^)=log(l + Qe-f). 



(39) 



This was used in [H] . However the proportion of the prey 
mass that is absorbed by the predator is not exactly the 
same in each feeding event. Variability arises for example 
from the difference in digestion between predator species 
and also from the difference in organic composition of 
prey organisms. In this paper we will allow variation by 
replacing the delta function by a Gaussian. So we will 
set 



a{y,z) = g a Az -ip(y)). 



(40) 



For the reproduction function b(y, z) we will use the 
product of Gaussians 



b(y,z) =b g a „(y-v)g a (z- /x). 



(41) 



This gives a mean offspringiparent mass ratio of e M and 
a mean mass ratio between parent before reproduction 
and parent after reproduction of e" . 

Finding the correct values for the parameters requires 
a close investigation of the data and is outside the scope 
of this paper. However, for the purpose of the plots, we 
have chosen parameters that appear at least reasonable 
from a biological point of view. For example, the pre- 
ferred predatonprey body mass ratio is believed to be 
around 10 2 or 10 3 [T3], so we have chosen /3 = 5.5. In or- 
der to estimate v, we need the average weight loss caused 
by reproduction processes. The average gonadosomatic 
index (ratio between the gonadal weight and the body 
weight) is actually measured for fish and is rather vari- 
able. It is found to be on average around 0.1 or 0.2 [52] . 
thus we have chosen v — 0.2 so that the average frac- 
tion of weight loss due to reproduction is around 20%. 
We have set the logweight difference between offspring 
and parent /i to a small value around -8 or -9. For the 
standard deviations in the parameters we use 17,3 = 2.5, 
(7„ = 0.05 and ct m = 0.5, although a careful analysis of 
the data will be necessary to determine them properly. 
We have set the absorption efficiency to Q — 0.9 [33J 
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1.2 



1.4 



FIG. 2. Plot of the curves for uq and r given in ( 36 1 and ( 35 1 
with the choices (37}, (40| and glj, for £ = 0.25. The rest 
of the parameters are p 5.5, ap = 2.5, ay, = 0, it = —10, 
<7 M = 0.5, v = 0.2, (j v = 0.05 and d — 0. The region of allowed 
p is shaded. 



because respiration and other metabolic processes have 
been modelled separately. In previous work |18j . where 
these processes were not separated, the net absorption ef- 
ficiency was replaced by a conversion efficiency of around 
0.2. In most plots we will set the variability in Q equal 
to zero, as well as the mortality rate. 

In Figure[2]we have plotted the maintenance rate r and 
the steady-state density coefficient uq as functions of p for 
the above choices of the parameters and £ = 0.25. The 
allowed interval for p appears shaded in that figure. It is 
encouraging that the observed value p ~ 1 (TU1 H21 H3J is 
contained within the interval. 



VI. STABILITY OF THE STEADY STATE 

It has been observed via numerical simulations in 
[51 [T7] that the power-law steady state is not always sta- 
ble against small perturbations but rather that the sys- 
tem can undergo a bifurcation in which the steady state 
becomes unstable and a stable travelling wave solution 
emerges. This phenomenon was investigated analytically 
in [TH] through a linear stability analysis. We now per- 
form a similar analysis in our generalised model. 

The only other paper that we are aware of that in- 
vestigates the stability of the power-law steady state is 
[7] but it deals, for reasons of simplicity, with a model 
where the growth due to feeding is independent of the 
prey density, thus avoiding having to deal with the asso- 
ciated non-linear terms. 

In order to discuss stability analytically, we consider 
the particular case £ = 0. According to (20), this corre- 



sponds to a maintenance rate proportional to the weight 
and a mortality rate independent of the weight. This is 
not quite realistic, but simplifies the analysis consider- 
ably because it leads to translational invariance in both 
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FIG. 3. Same as Figure [2] but with £ = and for two different values of Q. All the remaining parameters remain unchanged 
except for p = —7.5. Shaded regions show the intervals where r is positive, which determines regions of allowed p. We find 
that both cp and cp — rp — d are positive as well within that interval, leading to stable values of Mo- 



time and logweight. This will allow us to use the stan- 
dard Fourier transform. 



A. Perturbation in tto 

Before we consider the general, weight-dependent per- 
turbation we take a look at a particular perturbation that 
affects only the overall population density uq. So instead 
of ( 29 ) we consider the solution 



t(>,t) = u (t)e- px 



(42) 



where we now allow uq to depend on time. Substituting 
this into the evolution equation ( 27 1 with £ = gives 



duo 
~dt 



= -c P ul + (cb - rp - d)u , 



(43) 



where cp and cb are given in eq. (31). The solutions to 



this differential equation depend crucially on the strength 
of reproduction relative to maintenance and mortality. If 
reproduction is weak, i.e., if cb < rp + d, then the non- 
zero fixed point (steady state) at u — (cb — rp — d)/cp 
is unstable. If, however, reproduction is strong enough 
so that 



cb > rp + d, 



(44) 



then the fixed point is stable. In between, the system 
undergoes a bifurcation at which there is a whole line of 
fixed points exactly when cb = rp + d. 

In particular, this observation shows that the model 
without reproduction can never be stable against a per- 
turbation in the overall population density u . This in- 
stability was already noticed in [T5] , where it was argued 
that it could be avoided by suitable modifications of the 
model at the ends of the size-spectrum (plankton dynam- 
ics, senescent death). The observation that reproduction 
can provide a stabilising effect is new to this paper. 




FIG. 4. Same as Figure [3] left, except that p = —10. Shaded 
region show the interval for allowed p but unstable no (note 
that both cp and cp, — rp~ d are negative within the interval). 



As we have seen, the stability against variations in uq 
requires that cp and cb — rp — d are positive. In Fig- 
ure [3] we plot these coefficients together with uo and r as 
functions of p, for two different values of Q. As discussed 
in Section [Vj the parameter p is only allowed to lie in a 
certain interval where the maintenance rate is positive, 
which appears shaded. Within this interval, both cp and 
cb — rp — d are seen to be positive, so steady state is 
stable in this case. 

The region where both cp and c-Q — rp — d are negative 
corresponds to an unstable steady state, as shown in Fig- 
ure [4] with the same parameters but p = —10. For the 
following plots we will choose p = 0.9 as a suitable value 
leading to a stable steady state (with p w —8) and for the 
unstable solution we will choose p = 1.1 (for p w —10). 



10 




FIG. 5. The combined spectrum A(fc). Remaining parameters 
are: ft = 5.5, ap = 2.5, Q = 0.9, a^, = 0, cr M = 0.5, v = 0.2 
and ov = 0.05. Note that p = 0.9 lies in the region of stable 
Mo, whereas the fixed point for uq is unstable when p — 1.1. 



Inset contains the region of small k. 



B. General perturbation 

We now consider more general, weight-dependent per- 
turbations. It is convenient to change to the new density 
v related to u through 



u(x,t)=u$e px v(x,t), 



(45) 



so that the steady-state solution is just a constant 
v(x,t) = 1. We add a small perturbation e(x,t) to the 
steady-state solution 



v(x, t) = 1 + e(x, t) 



(46) 



and linearise the equation (27 1. Since the equation is 



linear and translationally invariant, we can solve it by 
means of the standard Fourier transform. 

We can express any perturbation as a linear combina- 
tion of plane waves labelled by a wavenumber k, 

e k (x,t) = e^ kx+ut l (47) 

In terms of that wavenumber, we get the following non- 
local dispersion relation 



iui{k) = -u dy s(y) 



e ikz , e -ik(y+z) _ ^ 



+ dy dz b(y, z) \ e - pv {e lky - 1) + e {f - 1)z {l - e- y )(e~ ikz - 1) 



(48) 



+ ikr, 



where we have used the steady-state condition ( 30 ) to eliminate the mortality rate d 



The sign of X(k) — —lm(u)(k)) determines stability. If X(k) is positive then the amplitude of the plane wave (47) with 



wavenumber k grows exponentially with time, rendering the steady state unstable. We find that X(k) = Xp(k) + Ab(/c) 
where 



and 



M*> - «. / * *> (-(✓■ + 1) «(*) + / .) («(».) + ~M + •» - 1)) 



A B (fc) = / dy ( dz b(y, z) (e~ py (cos(ky) - 1) + e ( ^ 1)z (l - e- y )(cos{kz) - l)j 



(49) 



(50) 



Note that the maintenance rate parameter r and the 
death rate parameter d no longer appear in these ex- 
pressions. Because a parent always uses weight during 
spawning, y is positive wherever b(y, z) is nonzero, and 
hence we can see that Ae(fc) is negative for any nonzero 
k. This shows that reproduction always has a stabilising 
effect. 

In the remainder of the section we will discuss the 
consequences on the stability of the steady state for the 
choices (37), (40) and (41) for the reproduction and pre- 
dation functions. We show in Figure [5] the combined 
eigenvalue spectrum A(fc) = Ap(fc) + Aa(fe) for two differ- 
ent values of p, corresponding to both a stable (p = 0.9, 



p = —8) and unstable (p = 1.1, p = —10) fixed point 
uq. We find that the spectrum for p = 0.9 is everywhere 
negative, corresponding to a stable steady state. We can 
see in the inset that the spectrum for p = 1.1 is posi- 
tive for small wavenumbers, leading to an instability of 
the steady state against very long-wavelength perturba- 
tions. At higher k the spectrum is more negative, in- 
dicating stronger stability against short-wavelength per- 
turbations. 

In Figure [6] we show the contribution from predation 
Ap(fc) for two different values of j3 with Q — 0.9 and for 
two different values of Q for ft = 5.5. Increasing Q from 
the value of 0.2 used in [THj has a considerable stabilising 
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FIG. 7. Effect of the variability in Q. Here we plot the differ- 
ence Ap(fc) — Xp - 1 (k), where A^f^fc) is the value for = 0. In 
the inset we plot a zoom of the spectrum for low wavenum- 
bers. Note that the magnitude of the effect of the variability 
in Q is very small. Remaining parameters are: /3 = 5.5, 
= 2.5, Q = 0.9 and p = 0.9. Inset contains the zone of 
small wavenumbers. 



effect. Note that we are allowed to increase Q because in 
this work we have separated out the losses due to main- 
tenance processes. We can see as well that decreasing 
the preferred body size ratio between predator and prey 
has a stabilising effect. Nevertheless, for realistic values 
of the parameters, the contribution from predation alone 
is positive at some wavenumbers k, showing that repro- 
duction is required to achieve stability. 

In order to characterise the effect of variability in Q we 
have considered non-zero in Figure[7j We have to keep 
the standard deviation sufficiently small so that the 
probability that the absorption efficiency is above 100% 
is negligible. This leads us to impose that 



< (1 - Q)e 



(51) 



for some typical value y of the log of the predator:prey 
mass ratio. In practice, typical values are 1 — Q « 0.1 
and y sa 5, so <C 6 • 10 -4 . Therefore the fact that the 
predatonprey mass ratio is so large implies that has 
to be very small, and to check its influence in stability we 
have chosen values around 10 -4 . The effect is negligible 
for small wavenumbers, although become slightly appre- 
ciable for highly oscillating plane waves. In Figure [7] we 
plot the difference Ap — Ap " 1 , being \^ the real part of 
the eigenvalue for = 0. Although the effect is very 
small, the variability in the feeding efficiency always en- 
hances the stability of the steady state. In Appendix [D] 
we show that the effect of these small fluctuations in the 
feeding efficiency consist on adding a diffusion term to 
the model. 

We have also studied the contribution Ab that repro- 
duction makes to the eigenvalue spectrum for various val- 
ues for /j, and v. The results are shown in Figure [8] As 
explained earlier, the spectra are always negative, show- 
ing that reproduction has a stabilising effect. Asymptot- 
ically Ab converges to a negative constant. Variations in 
fi and v affect the oscillatory behaviour found for small 
values of k. 



VII. CONCLUSIONS 

In this paper we made use of the fact that the power- 
law size spectrum that is observed in the pelagic ecosys- 
tem will be predicted by any dynamic model that is in- 
variant under scale transformations. That allowed us to 
generalise earlier models without spoiling the prediction 
of a power-law steady state. Where earlier models only 
included the effects of predation and intrinsic mortality, 
we added terms modelling maintenance costs and repro- 
duction and also allowed variability in the absorption ef- 
ficiency. Inclusion of maintenance and reproduction has 
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increased the stability of the steady-state solution. 

We did not go into much ecological detail in this paper 
and made no attempt to determine the parameters of 
the model directly from ecological data. Nevertheless, by 
exploiting scale invariance, we made several observations 
that are of ecological relevance and that were not clearly 
made in previous work: 

1. The power-law exponent for the size spectrum is 
fixed solely by the scaling properties of the parameters of 
the model. No detailed investigation of the model and its 
solutions is required to determine it. The exponent does 
not depend on details like the preferred predator:prey 
mass ratio, the feeding efficiency, the variability in feed- 
ing behaviour, the absorption efficiency, the maintenance 
costs, the mortality rate, or the details of reproduction. 

This is in contrast to the results of earlier works in 
which only predation was considered. In that special case 
there are not enough scaling relations to fix the steady 
state exponent and it will depend on the details of the 
model. 

It is a crucial aspect of our model that it contains 
both processes that are density-dependent (predation) 
and processes that are density-independent (maintenance 
respiration, intrinsic mortality, reproduction). It is the 
relative scaling of the rates for these processes that de- 
termines the steady-state power-law exponent. 

Camacho and Sole [5j studied the steady-state power- 
law exponent in a model with intrinsic mortality and re- 
production. However they assumed that mortality and 
reproduction rates were proportional to predation rates. 
Thus, in effect, all their processes were assumed to be 
density-dependent and again the steady-state exponent 
was not determined by scaling arguments alone. 

2. The assumption of scale invariance leads to pre- 
dictions about the scaling behaviour (151 of the various 



parameters of the model, and these can be tested through 



We can also make a prediction that has not yet been 
tested. The prey selection function s can be related to 
data on the stomach contents of fish as follows. Let 
l(x,x')dx' be the observed average number of prey with 
logweight between weight x' and x' + dx' found in the 
stomach of a fish of logweight x. This stomach content 
reflects what prey a fish has been eating recently, which 
is determined by the same predation rates that we used 
in constructing our model. Thus we get (see equations 
((221) and ([23])) 



l(x, x') 



s(x — x')u(x')T, 



(52) 



where T is a time related to the speed of digestion. We 
can assume that the fish density u(x') is close to the 
steady-state density (29 1 . Substituting the steady-state 
density u(x') 



u e 



- px 



into the above equation gives 



l(x,x') =UQe-^ x s{x-x')T. 



(53) 



So scale invariance makes a prediction for the form of the 
stomach content data. There is a lot of data available [34] 
and it should be possible, via a careful analysis of this 
data, to determine the scaling exponent £. We predict 
that this will confirm that £ « 1/4. 

3. From the condition that in the steady state the 
number of individuals must be conserved we deduced a 
relation between the parameters of the model which in 
particular restricts the exponent p to an interval. We saw 
that for biologically reasonable values of the parameters 
this interval is around p ~ 1, which is in agreement with 
observations. 

4. We have seen that the steady state in the model 
without reproduction, and with £ = 0, is unstable against 
perturbations in the overall population density uq. In 
fact this was our motivation for including a reproduc- 
tion term in the model in this paper. We have derived 



observation, as discussed in Section III 



the condition ( 44 ) for the magnitude of the reproduction 
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rate with respect to the magnitude of the maintenance 
and mortality rates which ensures stability against this 
perturbation. 

5. We have studied the stability of the steady state 
against all small perturbations in the case £ = 0. The 
stability is determined by the sign of X(k), which has the 
two contributions Xp(k) from predation and Ab(&) from 
reproduction. Maintenance and mortality do not enter 
these expressions. The contribution Xp(k) coincides with 
that calculated in [TSJ in the case of fixed absorption effi- 
ciency, however with the conversion efficiency K replaced 
by the much larger absorption efficiency Q. We have seen 
that this has an important stabilising effect. The contri- 
bution Ab(&) from reproduction is always negative, thus 
enhancing stability. 

6. We have generalised previous models for the pre- 
dation process to allow for variability in the absorption 
efficiency and have found that this does not have a big im- 
pact on the stability of the steady state. In Appendix |D| 
we show that this small effect can also be modelled by an 
additional diffusion term in the equations, which leads to 
simpler expressions. 

In future we intend to extract further interesting infor- 
mation from the model by taking more detailed ecological 
information into account. In particular we intend to use 
theoretical results from metabolic theories about the en- 
ergy budget [22[35l[36] to determine the relative strength 
of the various processes in our model and observational 
results to choose appropriate values for parameters. 

The analytic derivation of solutions and study of their 
stability that we performed in this paper should be 
pushed further. Numerical studies of the jump-growth 
model [17] has hinted at the existence of travelling wave 
solutions. It would be nice to learn more about them 
within the framework of scale-invariance, possibly by ex- 
tending our analysis in Appendix [C] We would also like 
to investigate the speed at which perturbations move 
through the size spectrum. 

Our stability analysis was restricted to the case £ = 0. 
In this case the symmetry algebra Aff(l,K), generated 
by time-translations and scale transformations, simplifies 
to the abelian symmetry group of translations in t and 
x. This allowed us to perform the stability analysis in 
terms of plane waves. It would be nice if the technique 
could be extended to the non-abelian case, possibly using 
the techniques from non-commutative harmonic analysis 
[37]. 

The power-law size spectrum has been observed over 
many orders of magnitude and covers not only fish but 
also all types of plankton and even inanimate particles. 
Our model is appropriate only for organisms that feed 
by swallowing smaller organisms and that reproduces by 
spawning a large number of smaller organisms. It is 
not appropriate for phytoplankton or inanimate particles. 
Other models are needed for these and it will be interest- 
ing to see how these models can be coupled together. We 
propose that again the guiding principle should be scale 
invariance. 



There is a limit to the amount of detail that can be 
incorporated into a community size spectrum model like 
the one described in this paper. In particular, because 
only the size of individuals is taken into account, their 
different behaviour in different life stages can not be 
modelled. A first step in the direction of refining the 
model was taken in |38) where an individual is described 
not only by its size but also by one species-specific trait, 
namely its size at maturity. That allowed individuals to 
follow more realistic ontogenetic growth curves. An in- 
teresting observation of [38] was how the size spectra for 
the different species, each of which singles out a partic- 
ular scale in the form of the maturity size, combine into 
a power-law community spectrum described by a scale- 
invariant McKendrick-von Foerster equation. 

In [35] the size of offspring was encoded through a 
boundary condition on the species' size spectrum but 
when combining these spectra into a community size 
spectrum the boundary condition had to be ignored in 
order to achieve a power-law steady-state size spectrum. 
It would be interesting to see if that work could be ex- 
tended to include a realistic model of the reproduction 
process and still produce a power-law community spec- 
trum. 

In a pure size spectrum model like ours, that does not 
take life stages into account, a realistic modelling of re- 
production is not possible. For example, the offspring of 
many marine species start their life in an egg stage dur- 
ing which they are not part of the consumer spectrum 
but are already preyed upon. In our model we had to 



neglect this fact. The scaling relation (19 1 implies that 



the typical size of offspring is proportional to the size of 
the parent. However, because we can not model the egg 
life stage, we can not claim that this gives a reliable pre- 
diction for the relation between egg size and parent size. 
Indeed, while such a relation may hold for copepods [39] , 
it certainly does not hold for fish [40 , where most species 
lay eggs of a similar size. 

One might be tempted to simply replace the reproduc- 
tion function (191 by one that represents the fact that 



most fish lay eggs of a similar size. However this would 
immediately lead to a steady state solution that is not 
a power law. The steady state spectrum would exhibit 
a peak at the preferred egg size which is not observed 
in the size-spectrum data. The correct approach is not 
to break the scale invariance of the community spectrum 
model, but instead to study how more detailed models 
at the species level can combine into a scale-invariant 
model at the community level, in the spirit of [381 141) . 
The species-level models would have to take into account 
effects like the duration of the egg stage, the separation 
between spawning grounds and feeding grounds, the var- 
ious spawning strategies, like for example the laying of 
eggs in clusters, and many more. How and why the com- 
bination of many such detail-rich species-level models can 
lead to a scale-invariant community-level model is an in- 
triguing problem and it is to be hoped that the experi- 
ence that mathematical physics has with the emergence 
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of scale invariance in complex systems can be exploited. bution 11(77, t) — ^ 1 ^ 2 -P(n, t) and 
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Appendix A: Systematic expansion of the master 
equation 

Our model is a Markov model and can be described 
by a master equation that gives the time evolution of the 
probability -P(n, t) that the system is in the state n at 
time t. There will be a contribution from each of the 
processes involved, 



dP _ (dP 
~dt ~ (~dt 



dP 
~dt 



dP 
~dt 



dP 
~dt 



(Al) 



This appendix will be devoted to the systematic expan- 
sion of the master equation in powers of the inverse 
system volume f2 _1 [20] . This expansion separates the 
macroscopic behaviour from the fluctuations, and gives 
a linear Fokker-Planck equation describing the fluctua- 
tions. Since the procedure is quite similar for the four 
processes we are considering, we will get into details for 
the predation process and will give the contributions of 
the other processes thereafter. An alternative derivation 
of the deterministic, macroscopic equation can be found 
in0. 

A concise way of writing the contribution of predation 
events to the master equation uses the step operator no- 
tation [20] . Since the probability to undergo a predation 
event in dt is Pijknirijfl -1 , we have 

(rair)- = £ ^w^ 1 -mn^p^t)}, 



hj,k 



(A2) 

where the step operators act on any function /(n) as 
EJ(n) = /(..., m + 1, . . . ) and E~ 1 /(n) =/(..., n, - 
!)■■■)■ 

In Section [TT] we already introduced the split of each 
variable rii(t) into a deterministic, macroscopic compo- 
nent <f>i (t) describing the density of individuals in weight 
bracket i, and a fluctuation component rji(t) as 



(A3) 



dP(n,t) 



-l/2 dIL (V>t) 

dt 



^ dn(T7,t) 



drji dt ' 



(A4) 

where we have used that VL 1 ' 2 drj i /dt — —d^ijdt (this 
follows from ([I]), taking the time derivative for fixed rij). 
The step operator E, now transforms 77, to rji + 0~ 1//2 , 
and can be expanded as 



E; = I + £T 1/2 



d_ 
dm 



2 H 



(A5) 



Substituting this expansion into the master equation 
( A2 ) we arrive at an equation containing different powers 
of the system volume £7. The highest order (S7°) terms in 
the expansion contain only macroscopic variables 4>i & n d 
vanish if these satisfy the macroscopic equation 

Terms at next order (S! -1 / 2 ) give a linear Fokker- 
Planck equation for the probability distribution 11(77) °f 
the fluctuations, 



an 

~dt 



- Va 1 

9 *: 



d 2 n 



By introducing the symmetric combination 
1 



fijk — [Pi 



ijk 



P, 



jik) 



(A6) 



(A7) 



we can give concise expressions for the coefficients in the 
Fokker-Planck equation, 



N 



I 



flji<f>l) , 

- fjli<f>j<l>l) , 

- fuj<f>i<t>i - fijiMj) ■ 



(A8) 



Terms at higher order in f2 -1 specify how the fluctuations 
deviate from being Gaussian. Fluctuations are seen to be 
damped by a factor fi -1 / 2 , and the volume fl of our sys- 
tem is very large. This justifies that we can concentrate 
on the study of the macroscopic equation in this paper. 

The remaining processes involved in the model can be 
treated in a similar fashion. The master equation for 
reproduction can be written as 



( dP(n,t) 
\ dt 



B life (E i E7 1 E A 



The new stochastic variables rji have a probability distri- 



' -I)[n,P(n,()]. 

(A9) 

At leading order the expansion in powers of f2 gives the 
deterministic contribution Q. At next order the fluc- 
tuations are governed by a Fokker-Planck equation like 
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( A6 ) but with coefficients given by 

Jl 

Lfj = ^ + B jU m jli) > 

I 

jl 

N$ = 22 {-Bijifa - /»',;, o, + Ih ,,in i .,<:>,) 



The effect of maintenance on the time evolution of the 
probability P(n, t) is given by 



(A10) 



(dP(n,t) 
\ dt 



^RifE^jEi-^lmP^t)], (All) 



which leads to ([5| at leading order and to the follow- 
ing contributions to the coefficients of the Fokker-Planck 
equation, 



Lu — —Hi, 

Lfj = S i+ ijRj, 

N% = Ri<f>i + Ri + icp i+1 , 



(A12) 



Nf, 



Finally, the contribution of intrinsic mortality to the 
master equation is 

(^^) D =X)A(E i -I)[n i P(n,t)] J (A13) 

which gives ^ and for the fluctuations we get contri- 
butions to the diagonal coefficients in the Fokker-Planck 
equation, 



Lfi = -A 
NS = Didn 



(A14) 



and the contributions to the off-diagonal coefficients are 
equal to zero. Note that the master equation for the 



complete model ( Al I is simply the sum of all the contri- 



butions, so the combined effect of all the processes in the 
macroscopic and the Fokker-Planck equation is just the 
sum of all the terms. 



Appendix B: Positivity of uo 

We are going to show that under very natural assump- 
tions the steady-state population coefficient uq is positive 
when mortality is not too large and the steady-state ex- 
ponent satisfies p > 1/e — £. This inequality holds for 
the observed values that are in the region of p rs 1 and 
£ s» 1/4. The assumptions are that 

1. Predators only eat prey that are smaller than them- 
selves. This means that s(y) ^ only if y > 0; 

2. Predators always gain weight during feeding, i.e., 
a (y, z ) ^ only if z > y; 

3. Offspring are always smaller than their parent and 
the parent always looses weight during spawning, 
i.e., b(y, z) ^ only if both y > and z < 0; 

4. All the functions s, a and b are non-negative and 
have non- vanishing integrals. 

We define the functions q%(p) = cp + (p + £)/p and 
<h (p) = c b — d + (p + £)/b- These are the denominator 
and numerator in the expression (36) for uq. We will 



show that these two functions are both positive under 
the above assumptions, which implies the positivity of 
Mo- 

The derivative of qi with respect to p, keeping the other 
parameters constant, is 



ggi 
dp 



dys(y)ePy ( y + (p + £ - 1) / dz(y + z)e^ z a(y, z) ] + /, 



(Bl) 



Assumptions 1 and 2 imply that, wherever the integrand 
is nonzero, we have z > y > 0, hence e^ p+ ^ z > 1 and 



dz(y + z)e < - p+ ^ z a(y, z) > y, 



(B2) 



since a is normalised to unity. Substituting this back into 
(Bl I and using assumption 4 shows that dqi/dp > for 
all p. Moreover, at the particular point p = — £ we can 



calculate 



qi (-£) = / dy s(y)e-t« > 0. 



(B3) 



Therefore the monotonicity of q\ implies that q±(p) > 
for all p > — £. 

Similarly we calculate 



16 



dp 



= I dy I dz b(y, z) (l - ye-^ + ^ + e~*(l - e- y ){ze {p+i)z + 1)) . 



(B4) 



According to assumption 3, y > and z < wherever 
the integrand is nonzero, therefore 1 — e~ v > 0. For 
p > 1/e— £ we also have ye~( p+ & y < 1 and ze (p+ ^ > -1 
wherever the integrand is nonzero. This, together with 
assumption 4, implies that dq2/dp > and hence q2 is 
strictly increasing with p. This implies the positivity of 
q-i for all p > 1/e — £ if (72 is nonnegative at p = 1/e — £, 
i.e., if 

< C B | p= i/ e -£ + -/B|p=l/e-$- ( B5 ) 

This represents a positive bound on d because /b is al- 
ways positive and cb > at p+£ = 1/e. In the particular 
case of absence of intrinsic mortality (d = 0) the above 



restriction automatically holds. We have proven that uq 



is positive for p > 1/e — £ if (B5 1 is true 



Appendix C: General scale-invariant solution 

In this appendix we will look at solutions that are scale- 
invariant but not time-independent. They take the form 

u(x, t) = e- px f( X ) where x = x- ln(t) /£ (CI) 

for some function /. Substituting this Ansatz into the 
evolution equation (27) gives 



■je^f'(x)= I dys(y) - e ra /(x)/(x - v) ~ e~ iy f(x)fix + v) + / dz a(y, z)e^>*+'V(x - *)/(* ~ V - z) 



+ I dy j dzb(y,z)[~ f( X ) + e-^yf(x + y) + e^-P-^(l-e-y)f(x + z) 
(p + 0f(x) + f'(x)} -df(x)- 



(C2) 



Unfortunately, this ordinary integro-differential equation 
for f(x) is still difficult to solve in general. 

An analytic solution can be found in the special case 
where only predation is considered. It is given by 



Six) = /oe ?x . 
where the prefactor /o is determined by 



(C3) 



l//o= / dys(y)(e^y + l- I dza(y,z)e^^). 

(C4) 

In terms of u this solution reads 



u(x,t) = f e 



— f nP -(p-€) x 



It. 



(C5) 



Note how the exponents p and £ only appear in the com- 
bination p — p — £ whose value is determined by the 
scaling behaviour of the feeding function alone, see (15). 



The model without reproduction, maintenance and in- 
trinsic mortality was treated in earlier work [21 [6l [8] , but 
this time-dependent power-law solution is new. 



Appendix D: Variability in absorption efficiency 

In the existing literature [HI HE], the absorption effi- 
ciency was taken to be the same in all predation events. 



In Section V B we were allowing variability in the absorp- 
tion efficiency, modelled by a Gaussian distribution ( 40 ) 
with standard deviation ay,. The fact that the preda- 
tor:prey mass ratio is so large implies that has to be 
very small (a^ <C 6 • 10~ 4 ). This allows us to approxi- 
mate the integral containing the Gaussian by Laplace's 
method. We will see that this will lead to a diffusion term 
added to the model with fixed absorption efficiency. 
With a chosen as in ( 40 1 the third term in the feeding 



part ( 23 1 of the model reads 



F 3 = e<"-0* J dy s(y) J dz e~^>g^(z - i>(y)) 

u(x — z)u(x — y — z) . (Dl) 

After a shift in the integration variable z this takes the 
form 

' f dy ( dzh(x,y,z)e- z2/2,T ^, (D2) 



^3= ... 
where 

h(x,y,z) = e^-^ x - z -^s{y) 

u(x - z - i/)(y))u(x - y - z - ip(y)). (D3) 

We can expand h(x,y,z) in Taylor series of z around 
z = up to second order and use Laplace's method (see 
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e.g. [H]) to evaluate the asymptotic behaviour of the 
integral. Neglecting exponentially decaying terms, we 
can approximate 



F, 



dyh(x,y,0) + 



dy 



dz 2 



h(x,y,z) 



in the limit of a^p <C 1. Taking into account that 

d 2 d 2 
dz- 2f{x - z) = dx 2f{x - z) 



(D4) 



(D5) 



holds for any sufficiently smooth function /, wc finally 
get the asymptotic expansion of the model for u^,<l, 



E(x), (D6) 



where 



J(x) = e {p -t )x J dys(y) (u(x - y) + e {p -^ )v u(x + y) 

(D7) 

E{x) = e^-V* J dys(y)e- ( - p '^ < - y) u(x - tp(y)) 

u(x-y-ip(y)). (D8) 



Therefore, allowing small fluctuations in the feeding ef- 
ficiency has the effect of adding a diffusion term to the 
model. 
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